library(Biostrings)
library(XML)
library(RMySQL)
library(GenomicRanges)
library(ggbio)

baseDir <- "/home/haseong/alpha/dev/FT-NGS"
targetDir <- "/home/haseong/alpha/dev/FT-NGS"
setwd(baseDir)

file.names<-dir(pattern="^KHS.*fna$")
dir.names<-dir(pattern="^out_KHS.*fna$")

fn<-file.names[4]
contig.file<-paste(baseDir, "/", fn, sep="")
contig.seqs<-readDNAStringSet(contig.file, "fasta")
out.dir<-paste(targetDir, "/out_", fn, sep="")
anno.file.names<-dir(out.dir, pattern="^anno.*RData")

for(cs in contig.seqs){
  contig.name<-unlist(strsplit(names(cs), split=" "))[1]
  anno.file.name<-paste("anno_", contig.name, ".RData", sep="")
  cs.anno.file.name<-paste(contig.name, ".RData", sep="")
  cat("Loading annotation file...\n")
  ## read anno.matrix
  load(paste(out.dir, "/", anno.file.name, sep=""))
  ## read blast result (bl)
  load(paste(out.dir, "/", cs.anno.file.name, sep=""))
  
  ##bl[1,8:9]
  gi.numbers<-unlist(lapply(strsplit(as.character(bl[,3]), split="\\|"), function(x){return(x[2])}))
  gi.number.names<-names(table(gi.numbers))
  gr<-list()
  cnt<-1
  for(gi in gi.number.names){
    sel.anno.matrix<-anno.matrix[which(anno.matrix[,1]==gi),]
    if(nrow(sel.anno.matrix)>0){
      sel.anno.matrix.sort<-t(apply(sel.anno.matrix[,5:6], 1, function(x){sort(as.numeric(x))}))
      ex.idx<-which(is.na(sel.anno.matrix.sort[,1])|is.infinite(sel.anno.matrix.sort[,1])|is.na(sel.anno.matrix.sort[,2])|is.infinite(sel.anno.matrix.sort[,2]))
      sel.anno.matrix.sort<-sel.anno.matrix.sort[-ex.idx,]
      sel.anno.matrix<-sel.anno.matrix[-ex.idx,]
      ir<-IRanges(sel.anno.matrix.sort[,1], sel.anno.matrix.sort[,2])
      ## add meta data
      ##df <- DataFrame(annotations=sel.anno.matrix[,3])
      ##values(ir)<-df
      ##autoplot(ir)
      ##names(ir)<-sel.anno.matrix[,3]
      ##gr<-GRanges(seqnames=sel.anno.matrix[,3], ranges=ir)
      
      sel.bl<-bl[gi.numbers==gi,]
      sel.bl.sort<-t(apply(sel.bl[,8:9], 1, sort))
      qir<-IRanges(sel.bl.sort[,1], sel.bl.sort[,2])
      ##names(qir)<-sel.bl[,3]
      mtch<-findOverlaps(qir, ir)
      mtir<-ir[as.matrix(mtch)[,2]]
      ## selected meta data
      mtanno=sel.anno.matrix[,3][as.matrix(mtch)[,2]]
      ##values(mtir)<-mtdf
      ##autoplot(mtir)
      gr[[gi]]<-GRanges(seqnames=Rle(contig.name, length(mtir)), ranges=mtir, anno=mtanno)
      ##qir[as.matrix(mtch)[,1]]
      ##subsetByOverlaps(ir, qir)
    }
    cnt<-cnt+1
    cat(gi, "/", cnt, "/", length(gi.number.names), "\n");flush.console()
  }
  
  grl<-GRangesList(gr)
  autoplot(gr[1:5], geom = "alignment", group.selfish = FALSE)
    
  
  
}





##exe.url<-"http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?"
##exe.query<-paste(exe.url, "db=nuccore&id=", gis[1], sep="")
##tmpxml <- getURL(exe.query)
}


bl<-list()
cnt<-1




library(NCBI2R)



#### =======
library(RCurl)
library(XML)
library(snippets)
library(tm)

# esearch
url <- "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?"
##q   <- "db=pubmed&term=saunders+nf[au]&usehistory=y"
q   <- "db=nucleotide&term=CP002177[accn]&usehistory=y"
##q   <- "db=nucleotide&term=325121063[gi]&usehistory=y"
esearch <- xmlTreeParse(getURL(paste(url, q, sep="")), useInternal = T)
ids  <- xmlValue(getNodeSet(esearch, "//IdList")[[1]])
webenv  <- xmlValue(getNodeSet(esearch, "//WebEnv")[[1]])
key     <- xmlValue(getNodeSet(esearch, "//QueryKey")[[1]])



# efetch
url <- "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?"
##q   <- "db=pubmed&retmode=xml&rettype=abstract"
q   <- "db=nuccore&id=325121063"
tmpxml <- getURL(paste(url, q, sep=""))
##tmpxml<-getURL("http://www.ncbi.nlm.nih.gov/nuccore/325121063")
##tmpxml <- getURL(paste(url, q, "&id=", ids, sep=""))
tmp<-xmlTreeParse(tmpxml, useInternal = T)
efetch <- xmlTreeParse(getURL(paste(url, q, "&id=", ids, "&WebEnv=", webenv, "&query_key=", key, sep="")), useInternal = T)
feat <- getNodeSet(tmp, "//FEATURES")

# words
abstracts <- sapply(abstracts, function(x) { xmlValue(x) } )
words <- tolower(unlist(lapply(abstracts, function(x) strsplit(x, " "))))

# remove parentheses, comma, [semi-]colon, period, quotation marks
words <- words[-grep("[\\)\\(,;:\\.\\'\"]", words)]
words <- words[-grep("^\\d+$", words)]
words <- words[!words %in% stopwords()]
wt <- table(words)
wt <- wt[wt > 1]
cloud(wt, col = col.br(wt, fit=TRUE))
